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1.  INTRODUCTION 


This  report  describes  progress  made  in  during  the  second  year  of  study  on  the  project  entitled  “Model  &  Expansion 
Based  Methods  of  Detection  of  Small  Masses  in  Radiographs  of  Dense  Breasts.”  Our  end  goal  is  to  develop  technology 
for  computer  aided  diagnosis  that  can  detect  masses  in  dense  mammograms  having  a  diameter  less  than  1  cm.  The  basic 
“idea”  of  this  IDEA  category  award  is  to  detect  subtle  masses  by  tuning  the  central  frequency  and  width  of  a  basis 
function.  By  modeling  the  shape  of  a  mass  through  this  flexibility  (i.e.,  changing  the  shape  of  the  “bump”)  we  hope  to 
better  detect  small  and  subtle  masses  in  dense  breasts  and  improve  the  chances  of  early  detection  through  screening 
mammography. 

This  annual  report  of  our  second  year  continues  progress  made  towards  accomplishing  the  goals  described  in  Tasks  la 
and  lb,  “Feasibility  assessment  and  design  of  model  based  method.”  In  addition  we  have  begun  work  towards  building  a 
detector  as  described  in  in  Task  2b,  “Detector  analysis  and  implementation.” 

During  the  past  year  we  have  focused  on  an  interesting  family  of  functions  that  are  well  suited  to  the  problems  of  mass 
detection:  fractional  B-splines  are  an  extension  of  the  B-spline  functions.  In  this  report,  we  study  the  influence  of  the 
spline  parameter  in  the  analysis  of  mammograms  and  we  describe  experiments  carried  with  aim  of  identifying 
parameters  most  efficient  in  characterizing  the  masses  in  dense  breasts. 

Recall  that  in  our  previous  report,  we  observed  that  masses  within  mammograms  are  represented  low  frequency  signals.. 
Thus  the  analyzing  functions  we  use  should  to  be  smooth  in  terms  of  shape  and  can,  for  instance  be  first  or  second 
derivative  functions.  The  fractional  spline  family  is  exactly  what  we  want  since  it  has  one  continuously  varying  order 
parameter.  Thus,  one  can  tune  the  central  frequency  of  a  wavelet  that  will  be  efficient  in  the  detection  of  masses  of 
arbitrary  size. 


2.  BODY 


It  has  been  shown  that  the  dyadic  transform  is  not  sufficient  for  the  case  of  mass  detection  [1].  Thus  we  have  constructed 
a  2D-continuous  discrete  transform,  using  the  “A  trou  Algorithm”  in  which  we  incorporated  the  fractional  spline 
wavelet.  In  addition,  we  will  report  on  a  component  of  a  detector  algorithm  which  uses  a  dictionary  of  bases  computed 
from  Wavelet  Packets.  Our  goal  related  to  Task  2,  was  to  find  an  “average”  best  basis  according  to  minimizing  a 
criterion  of  energy  within  a  wavelet  packet  expansion. 

In  this  report,  we  have  processed  numerous  Regions  of  Interest  (ROI)  to  allow  us  to  compute  this  average  best  basis  and 
extract  some  specific  features  for  the  detection  of  small  masses. 

Design  and  Methods 


The  2D  fractional  B  spline  Wavelet  Transform 
The  B  spline  functions 

Splines  [2]  are  piecewise  polynomials  that  are  smoothly  connected  together.  The  joining  points  of  the  polynomials  are 
called  knots.  For  a  spline  of  degree  n,  each  segment  is  a  polynomial  of  degree  n,  which  suggests  that  we  need  (n+1) 
coefficients  to  describe  each  piece.  However,  there  is  additional  smoothness  constrain  that  imposes  continuity  on  the 
spline  and  its  derivatives  up  to  order  (n-1)  at  the  knots,  so  that,  effectively,  there  is  one  degree  of  freedom  per  segment. 
Here,  we  will  only  consider  splines  with  uniform  knots  and  units  spacing.  The  remarkable  result  is  that  these  functions 
are  uniquely  characterized  in  terms  of  B-spline  expansions: 

six)  =  '^cik)/3"(x-k)  (0.1) 

teZ 
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which  involve  integer  shifts  of  the  B-spline  of  degree  n  denoted  by  ;  the  parameters  are  the  B-spline  coefficients 
c(k).  B-splines,  defined  below  are  bell-shaped  and  constructed  from  the  (n+1)  fold  convolution  of  a  rectangular  pulse 
(centered): 


y9°(x)  = 


1, 

]_ 

2’ 

0, 


1  1 

—  <  ;c<  — 

2  2 

otherwise 


(0.2) 


The  fractional  B-spline  functions  [2,3] 

We  can  extend  B-Spline  to  all  fractional  degrees  OC  >  •  These  splines  are  constructed  using  linear  combinations  of 

(X  I 

the  integer  shifts  of  the  power  functions  (one-sided)  or  pel  (symmetric): 


ur  = 


;c>0 

0,  otherwise 

\o: 


-2  sin 


^  n:  ^ 

—a 

2 


X^"  log  X 


/  ^  \l+n  ’ 

.(-0  ^ 

In  each  case  they  are  a-Holder  continuous  ford^  >  0 .  They  satisfy  most  of  the  properties  of  the  traditional  B-splines,  as 
described  below.  In  particular  the  Riesz  basis  condition  and  the  2-scale  relation,  which  make  them  suitable  for  the 
construction  of  a  new  families  of  wavelet  bases. 


a  =  2n  +  l 


a  =  2n 


Definition  and  notation 

Euler’s  gamma  function  is  defined  as: 

r(M  H- 1)  =  f  x''e~''dx  (0.3) 

Jo 

In  particular,  r(n  -h  1)  =  n!.  This  suggests  the  following  generalization  of  the  Binomial  coefficients 

T{u  +  \) 


In  particular,  this  definition  implies  that 


(0.4) 

r(v-i-i)r(M-v-f-i) 

=  0  when  k  is  a  strictly  negative  integer.  Moreover,  for  W  >  0 ,  we  have 


the  well-known  binomial  theorem:  (l  + 


K*. 


The  differentiation  operator  can  be  extended  to  non-integer  exponents  rather  simply  in  the  Fourier  Domain; 
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D'' f  {x)  <r^  {j(Of  f  {(O)  where  fiCO)  =  I  f(x)e  denotes  the  Fourier  Transform  of  f  (x)  and  where 

j  —  ylHi  arg(z)G  ^].  This  is  equivalent  to  Liouville’s  definition  of  fractional 

derivative  operator. 

The  forward  fractional  finite  difference  operator  of  order  a  is  defined  as 

A:f(x)=f^i-iy(^]f(x-k).  (0.5) 


k^O 


It  is  a  convolution  operator  which  has  a  more  straightforward  interpretation  in  the  Fourier  domain 


A:(fy)  =  (l-e-^")“=5;(-l)* 


^=0 


\h 


-jcok 


(0.6) 


The  natural  building  blocks  for  the  fractional  splines  are  the  one-sided  power  functions  x^ ,  which  have  precisely  one 
singularity  of  order  a  (Holder  exponent)  at  the  origin.  Their  Fourier  transform  is  r{(X  +  l)/(jCOy^^ . 

Similar  formulae  exist  for  the  symmetric  functions  whose  Fourier  transform  is  r(61f  H-l)/|7^y|  .  However,  in  our 

study,  we  only  consider  the  one-sided  functions. 


Fractional  B-splines 

We  consider  here  the  causal  fractional  B-splines.  By  analogy  with  the  classical  B-spline  we  define  the  fractional  B-spline 
by  taking  the  {cc-\-\^th  fractional  difference  of  the  one-sided  power  function 


Aflf  a  i  +00 

p-  (X) = = — - —  y  (-1)* 

r(a+l)  na+Dt'o 

1 


f  a+l^ 


(x-kt. 


(0.7) 


These  functions  are  for  61^  >  ~— .  While  they  seem  to  be  decaying  reasonably  rapidly,  they  are  not  compactly 

supported  unless  a  is  an  integer,  in  which  case  we  recover  the  classical  B-splines.  In  general,  they  do  not  have  an  axis  of 
symmetry  either.  The  fractional  causal  B-splines  satisfy  the  convolution  property 

(0.8) 


The  properties  of  fractional  splines 

The  Fourier  transform  of  the  causal  fractional  B-splines  are  given  by 

_.v..N«+l 

p:m= 


V 


JOi 


(0.9) 


This  equation  is  indeed  compatible  with  the  convolution  property. 

Fractional  derivatives  and  regularity 

One  of  the  primary  reasons  for  the  success  of  polynomial  splines  in  image  and  signal  processing  applications  is  that  they 
can  be  differentiated  very  simply  by  taking  finite  differences.  This  property  generalizes  nicely  to  the  fractional  case: 

D^p:{x)=s\pr{x)  (0.10) 

(The  2  scale  relation  is  described  in  the  next  section). 

Fractional  spline  spaces  and  Riesz  bounds 

The  mathematical  space  of  fractional  splines  of  degree  a  with  knots  at  the  integers  is  defined  as: 
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=  (x-k):cElA. 

[  keZ  ) 


Proposition:  For  00-1/2  the  fractional  B  spline  generates  a  Riesz  basis  of  .  Specifically,  one  has  the  following 
/2-L2  norm  equivalence 


Y,c{k)P'^  (x-k) 


Where  \  > 


(2\ 


yn  j 


keZ 

2a+2 


and  B^<l  +  2C(2a+2) 


71 


This  result  ensures  that  the  B-spline  representation  is  stable  and  that  fractional  spline  spaces  are  well-defined  (closed) 
subspaces  of  L2.  Finally,  the  fractional  splines  of  degree  a  have  a  fractional  order  of  approximation  a+1. 


Detinition  of  the  fractional  spline  wavelet 


Let  us  now  show  explicitly  how  to  use  these  B-splines  to  obtain  new  wavelet  families  with  a  continuously  varying  order 
parameter.  We  will  concentrate  on  the  case  of  orthogonal  wavelet  bases.  These  are  obtained  by  orthogonalization  of  the 
fractional  B-splines.  The  key  quantity  is  the  fractional  B-spline  autocorrelation  sequence: 

a,(k)  :=  {P:(x),fi:(x-k))  =  Pl“^\k)  (0.11) 

Where  the  right  hand  side  expression  is  a  direct  consequence  of  the  convolution  relations.  Moving  to  the  Fourier  domain 
we  get 

kez'  keZ 

since  P“  {x)P^  {x  +  n)dx. 

neZ 


To  compute  A^{e^^)  we  choose  to  perform  an  infinite  summation  in  the  Fourier  domain  because  of  the  relatively 
convenient  expression  of  {CO)  .  And  we  express  the  orthogonal  scaling  function  as 


^{x)  =  ^{a^frpt{x-k) 


keZ 


The  corresponding  two-scale  relation  is 

^{xl2)  =  yl2Y,K  {k)^ix  -  k) 


(0.13) 


(0.14) 


keZ 


and  the  refinement  filter  is  given  by 


a+\ 


(0.15) 


The  corresponding  orthogonal  wavelet  filter  can  then  be  obtained  using  Mallaf  s  method: 

Gq(z)  =  z.H„(-z~^).  (0.16) 

An  overcomplete  representation  for  signals:  the  “A  trous  Algorithme”[3] 


We  have  shown  that  the  fast  dyadic  transform  was  not  sufficient  in  the  case  of  small  mass  detection  [1].  However,  it  is 
possible  to  sample  the  time-scale  space  in  other  ways.  The  Continuous  Discrete  Wavelet  Transform  (CDWT)  keeps  the 
dyadic  sampling  of  scale  (Discrete)  but  provides  a  quasi-continuous  sampling  in  time  (Continuous).  It  is  similar  to  a  fast 
biorthogonal  wavelet  transform,  without  downsampling  in  the  time  domain.  In  terms  of  representation,  this  transform 
falls  between  the  large  redundancy  of  the  CWT  and  the  strictly  non-redundant  DWT. 
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cd^{j,k)  =  CWT^ia  =  2\t  =  k) 
cd,U,2'‘)  =  DWT,ij,k) 

Proposition:  For  any  J  ^  0 , 

^^,+1  m=<3,  r«i  1/  - 

where  a.  and  d-  are  the  approximation  and  the  detail  coefficients  of  the  signal 

hj  =t  2^'-*  [h,  ]  and  g.  [g,  ]  (0.17) 

1\  and  g|  are  respectively  the  low-pass  and  high  pass  filters  linked  to  the  wavelet. 


X 


COj 


cd^ 


cd2 


cd^ 


Figure  1.  Overview  of  processing  to  compute  the  ‘a  trou’  algorithm. 


The  2D  Continuous  Discrete  Fractional  Spline  Wavelet  Transform 


Implementation 


As  described  above,  we  used  the  “a  trou  Algorithm”  to  perform  this  transform.  We  appled  the  low-pass  and  high-pass 
filters  respectively  on  the  rows  and  the  columns  as  illustrated  in  Figure  1.  We  worked  in  the  Fourier  domain:  CGq  is  the 

Fourier  Transform  of  the  image,  the  initialization  is  therefore  ensured  by  the  Fourier  operation.  After  the  filtering 
operations,  an  inverse  Fourier  transform  is  performed  to  obtain  the  coefficients  in  the  spatial  domain. 

Filters  1\  and  gjare  respeetively  linked  to  G^of  (1.15)  and  (1.16)  with  Go(e^‘“)  =  e” (-c“^®)  .  In 


H.  G 

order  to  have  normalized  Quadratic  Mirror  Filters  (QMF),  we  consider  and  , 

v2  v2 
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ROWS 


COL 


ROWS 


COL 


Figure  2:  Scheme  of  the  overcomplete  algorithm,  ;  Approximation  in  the  wavelet  basis,  c?-  :  Horizontal  detail 
in  the  wavelet  basis,  df :  Vertical  detail  in  tbe  wavelet  basis,  df  Diagonal  detail  in  the  wavelet  basis. 


Practical  considerations: 

For  the  computation  [0.6]  of  the  Fourier  transforms  of  the  fractional  B-splines  filters  we  need  to  evaluate  the 

k 

autocorrelation  filter  at  V  =  — ,  for  ^  —  1.  For  this  we  use  the  Poisson-equivalent  expression  of  (6)  i.e., 

n 

^1)3^  {(0  +  27tk)  .  More  precisely,  since  we  must  limit  this  summation  to  a  finite  number  of  terms,  we  evaluate 

iteZ 

using  the  following  asymptotic  equivalent: 


.  .  .|2a+2  - 

2a±l 


1  («+1)[5  +  2’'']  (a+l)(2«+3)l 


.  (0.18) 


2a+l  Tvr2or+2 


N 


An  accurate  analysis  of  this  expression  shows  indeed  that  the  remainder  difference  between  the  Ihs  and  the  rhs  is 

(  1  "i 

0  — 2^^  ,  which  ensures  an  accuracy  that  is  larger  than  200  dB  with  N=100  computed  terms,  and  for  all  values 

yN  j 

of  a  >  ~  0.5.  (The  code  implementing  this  is  provided  in  the  Appendix  of  this  report.) 


Moreover,  we  use  biased  coefficients  to  ensure  positive  coefficients,  and  display  a  gray  scale. 
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fractionnal  spline  filtersiresponse  in  frequency 


1.4 

1.2 


H:  low-pass 
_  G  high-pass 

_  abs(H)^+abs(G)^ 


140 


Figure  3:  Filter  design  as  a  function  of  alpha,  for  a  128  points  signal. 

The  filters  are  a  quadratic  mirror  filter,  which  are  very  desirable  for  signal  processing  -  it  allows  us  to  use  existing  tools 
that  can  be  applied  to  such  filters.  This  is  also  the  case  for  the  wavelet  packet  representation  described  below,  next. 

Best  basis  Search  for  Detection 

The  ID  wavelet  packet  representation 

Wavelet  packets  generalize  the  link  between  multiresolution  approximations  and  wavelets.  Previously,  we  noted  [5]  that 
a  space  of  a  multiresolution  approximation  is  decomposed  in  a  lower  resolution  space  and  a  detail  space 

by  dividing  the  orthogonal  basis  of  into  two  new  orthogonal  basis,  one  for  and  the  other  for  . 

But  here,  instead  of  dividing  only  the  approximation  space  to  construct  the  detail  space  and  a  new 

approximation  space  ,  we  want  also  to  divide  to  derive  new  bases!  The  following  theorem  allow  us  to  perform 
such  divisions: 

Theorem  :  ( Coiftnan,  Meyer,  Wickerhauser) 

Let  ^ —  2^ li)^  ^  be  an  orthonormal  basis  of  a  space  U j .  Let  h  and  g  be  a  pair  of  conjugate  mirror  filters. 
Define: 

e%{t)=^h{n\eft-2^n)  and  d],ft)  g[n\d ft -2^  n) ,  the  family 

n=-oo  n--o° 

is  an  orthonormal  basis  of  U  j . 
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Suppose  our  signal  is  approximated  at  the  scale  2^.  The  corresponding  approximation  space  is  V^and  admits  the 
following  orthogonal  basis:  (^  —  2^ n)|  ^ . 

We  can  easily  represent  the  recursive  splitting  of  vector  spaces  in  a  binary  tree  -  is  the  root  of  the  tree,  and  to  each 
node  (/,p)  we  associate  a  detail  space  Wf  that  admits  the  following  orthonormal  basis:  —  .At  the  root, 

we  have:  .  For  each  node,y-L  is  its  depth  in  the  tree  and  p  is  the  number  of  nodes  that  are  on  its  left  at  the  same 

depth.  Suppose  now  that  we  have  already  constructed  Wf .  The  construction  of  its  children  Wf  and  Wf  is  given  by 
the  splitting  relation  of  the  last  theorem: 

¥%(t)=^h[n]wf(t-2Jn) 


And  since  (t  —  2-'^'  n) ,  (t  —  2^^*  n)}  ^  is  and  orthonormal  basis  of  Wjf’  we  have: 


Following  this  recipe,  we  constructed  a  binary  tree  of  wavelet  packet  bases.  We  may  choose  to  divide  each  node  further 
(or  not)  so  that  we  obtain  a  binary  tree  where  each  node  has  either  0  or  2  children.  A  tree  is  called  admissible  if 

at  the  end  of  the  exspansion:  =  ©[  where  { j. ,  p.  are  the  leaves  of  our  tree. 

The  union  of  the  corresponding  wavelet  packet  bases  w  “2^' defines  an  orthogonal  basis  of 


=  :  it  is  called  a  wavelet  packet  basis. 


Figure  4:  Example  of  an  admissible  tree  computed  from  a  wavelet  packet  expansion. 
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Best  basis  search 


We  implemented  the  wavelet  packet  decomposition  in  two  dimensions.  We  can  imagine  Figure  4  above  with  4  children 
per  node  (one  approximation  node  and  three  details  nodes  for  each  expansion).  The  search  for  a  best  basis  uses  an 
additive  criterion  described  below. 


The  criterion: 

As  an  early  choice  in  our  experiments,  we  selected  the  following  Entropy  additive  criteria: 


”  X 

Ent  =  -'YjP^o%^{p,)  with  p,  =  — =(x,.^.) 

k=\  X.  ^ 


,  (0.20) 


where  X  is  the  image  matrix. 

Algorithm  in  2D: 

Initialization:  =  ¥j’JL’{P^k)e  -l] 

Por  7  =  -^max-l’  0: 

^  ^  J y/f  if  )  <  Ent  )  +  Ent  ( )  +  Ent  ( )  +  Ent  ( ) 

'  I otherwise 
p  =  0,...,V-\ 

The  selected  leaves  of  the  tree  are  the  best  basis. 


The  criteria  is  an  entropy-based  metric  that  has  been  used  in  application  of  image  compression:  Its  purpose  is  to 
concentrate  in  a  few  coefficients  most  of  the  information  of  an  image.  In  our  case,  the  highest  value  coefficients  can 
describe  the  mass  if  the  wavelets  are  well  adapted  for  detection  (Please  see  discussion  section  below). 

Implementation  notes 

We  needed  to  adapt  the  algorithm  of  the  decomposition  of  wavelet  packets  for  the  end  goal  of  detection.  The  2D  discrete 
decomposition  already  existed  in  a  Wavelab  library.  We  have  modified  this  available  code  for  our  purpose,  as  described 
in  the  Appendix  of  this  report  (cf  best  basis  directory).  In  particular,  as  we  did  not  want  to  downsample  coefficients  we 
needed  to  preserve  normalization  of  the  entropy  when  we  compared  the  entropies  of  the  father  to  the  children  nodes  in 
the  tree.  Thus  we  divided  the  children’s  entropy  by  a  factor  of  4. 

We  also  stored  the  different  coefficients  of  each  best  basis  for  reconstruction  and  visualization  purposes. 

Thus,  we  forced  the  first  node  to  the  value  1  in  order  to  allow  the  program  sb_Plot2dpartition.m  to  decompose  the  nodes. 
This  code  also  provided  in  the  Appendix,  under  the  heading  of  Best  basis  directory. 


Results 


Computation  of  the  2D  Wavelet  transform 

We  computed  along  the  dyadic  scales,  the  CDWT2D  for  different  parameter  values  of  a  (please  see  the  subroutine 
"'simu_alpha_scale.m'"  in  the  Appendix  for  details).  The  results  in  terms  of  detection  efficiency,  are  discussed  in  the  next 
section.  Pirst,  we  choose  a  range  of  a  values  that  made  sense  in  terms  of  detecting  masses  less  than  1cm  in  diameter. 
First  we  computed  an  expansion  for  a  phantom  image  of  a  mass  (Figure  6)  with  a  large  amount  of  Gaussian  noise  added 
(Figure  7).  This  experimental  phantom  was  used  to  validate  (debug)  our  algorithm  and  confirm  function. 


12 


alpha 


Figure  9:  Expansion  of  approximation  coefficients  of  the  CDWT2D  for  a  phantom  containing  a  noisy  mass. 


fractionnal  spline filters:response  infrequency 


Figure  10:  Filter  frequency  response  for  corresponding  values  of  selected  a. 

We  tested  the  transform  for  the  following  a  values  shown  in  Figures  10  and  11.  This  range  in  alpha  values 
was  selected  because  the  filters  had  smooth  slopes,  desirable  for  matching  the  shape  of  masses. 


fractionnal  spline  filters:response  in  frequency 


Figure  12 
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Approximation 


Vertical  detail 


Figure  16:  CDWT2D  expansion  of  case  02. 


Figure  18;  CDWT2D  expansion  of  case  04, 


Best  basis  search 


We  applied  the  best  basis  search  strategy  described  above,  on  images  extracted  from  the  database  of 
mammograms  available  on  the  University  of  South  Florida  website 

http://marathon.csee.usf.edu/Mammographv/Database.html.  We  wrote  a  routine  that  allows  us  to  choose 
regions  of  interest  within  these  mammograms  (512x512  pixels).  Each  case  contained  a  cancerous  masses 
and  their  features  are  detailed  in  the  following  table: 


Cases 

Properties 

02,04,05,06,08,14,17, 

18,21,25,26,27,28,29,30 

Type:  cancer 

Density:  4 

Lesion  type:  mass 

Mass  options:  Irregular  and  lobular 

The  returned  results  are  trees  that  correspond  to  the  best  basis.  For  practical  visualization  and  practicality, 
these  results  are  displayed  in  terms  of  2D  trees  that  are  computed  as  if  the  images  were  not  down-sampled. 


Low-high 

Low-  low 

High-high 

hllh 

hill 

hlhh 

hlhl 

Figure  19.  Display  format  of  filtered  coefficients  using  a  wavelet  packet  expansion. 

All  processed  images  were  512x512  matrix  size  and  the  computation  was  time-consuming  in  terms  of 
processing.  We  therefore  decide  to  limit  our  expansion  to  a  depth  of  4,  which  corresponds  to  a  scale  2'^=32. 
We  also  used  the  filters  shown  in  Figure  20  below  for  computation. 

fractionnal  spline  filters:response  in  frequency 


Figure  20,  Selected  filters. 
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Discussion 


Influence  of  the  spline  parameter 

For  ai<a2,  the  bandwidth  of  Hai  is  larger  than  the  bandwidth  of  Ha2  (please  see  Figure  3).  The  frequency 
resolution  is  therefore  less  for  Hcti  than  for  Hct2,  but  its  spatial  resolution  is  better  as  observed  visually  in 
Figure  8. 

We  also  note  some  shifting  effects  perhaps  due  to  using  a  non-zero  phase  filter  -  the  position  of  the  mass  is 
changing  for  higher  values  of  a  (as  seen  for  a>3  in  Figure  9). 

The  use  of  fractional  splines  is  well  justified  when  examining  Figure  9  -  one  can  note  the  contrast 
improvement  between  the  mass  and  the  background,  leading  to  a  simple  detection  task  of  the  mass. 

For  different  scales  and  a,  it  is  interesting  to  see  different  details.  For  instance,  for  the  image  shown  in 
Figure  14,  the  couple  (scale=4,  a=-0.1)  seems  to  extract  the  most  visual  information  for  the  mass  whereas 
the  couple  (scale=2,  a=-0.5)  gives  a  poor  result  in  terms  of  detection. 

A  useful  range  of  a  concerning  the  images  seems  to  be  within  the  interval  [-0.4,1].  And  among  these 
adjustable  values  of  a,  it  would  be  interesting  to  be  able  to  choose  automatically  the  best  one  that  gives 
high  contrast  between  the  mass  and  the  background  tissues  of  the  breast.  We  did  not  observe  a  single 
couple  (scale,a)  for  all  our  sample  cases.  However  we  identified  some  significant  trends.  We  see  in  Figures 
8  through  18  that  making  the  value  a  vary  within  the  right  interval  is  very  interesting  since  it  allows  us  to 
see  details  present  within  the  mass. 
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For  higher  values  for  a  we  observed  that  the  details  and  the  approximations  did  not  result  in  a  superior 
representation  of  the  mass  shown  in  Figure  10.  Beyond  these  values  the  results  are  not  interesting  in  terms 
of  detection  since  the  corresponding  coefficients  emphasized  the  texture  within  the  mass. 

Indeed,  with  a  continuous  analysis  framework,  which  would  allow  decomposition  on  voices  between 
scales,  we  may  obtain  a  richer  parameter  space  to  identify  the  best  basis  for  mass  detection.  We  discuss  the 
application  of  a  continuous  transform  at  the  end  of  this  section. 

In  the  search  for  the  best  basis,  let  us  first  show  the  link  between  the  criterion  metric  and  detection  goal. 
The  selected  coefficients  within  a  best  basis  minimize  the  entropy  of  the  image.  However,  in  our  case  it  is  a 
these  same  coefficients  that  mostly  represent  the  presence  of  a  mass.  Indeed,  the  wavelet  based  on  the 
fractional  splines  simplifies  the  mass  detection  problem  when  a  good  choice  of  a  is  made.  The  minimum 
entropy  means  the  maximum  amount  of  information  is  contain  within  a  minimum  number  of  coefficients. 
The  coefficients  returned  by  identifying  a  best  basis  represent  this  information.  As  mention  in  our  previous 
report,  mammograms  are  largely  low  frequency  images.  Thus  entropy  is  linked  to  the  shape  of  a  mass  since 
minimizing  the  entropy  leads  to  optimal  correlation  between  basis  coefficients  and  the  mass  signal  (as 
shown  in  our  previous  year’s  report).  To  demonstrate  this  point,  we  show  in  Figure  21  that  the  entropy  is 
lower  for  a  signal  with  a  mass  than  for  an  image  containing  high-frequency  features. 


Entropy=-6.6008e+009  Entropy=-9.9922e+008 

Figure  21.  Samples  images  of  different  entropy. 


We  performed  experiments  on  three  different  values  of  the  a  parameter  in  fractional  spline  wavelets.  The 
aim  here  was  not  to  perform  the  transform  for  every  case  but  to  see  the  similarities  and  differences  in  the 
selected  a’s,  and  to  obtain  a  range  for  computing  best  bases.  In  the  next  phase  of  the  project,  these  bases 
will  lead  to  a  classification  scheme  for  masses.  At  this  point,  we  have  convinced  ourselves  by  judging  the 
basis  in  terms  of  a  visualization  by  reconstructing  the  mass  from  its  best  basis. 

We  observe  that  the  expansion  is  often  refined  in  the  top-right  quadrant,  which  corresponds  to  a  first 
decomposition  of  low  frequency  information.  In  addition,  we  note  that  for  higher  a  values,  more  nodes 
were  expanded.  Thus  we  needed  to  expand  within  coarser  scales  to  minimize  entropy. 
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Minimum  decomposition 
a=-0.4 

Average  number  of  decomposed 
nodes:  8.2 


Minimum  decomposition 
a=0 

Average  number  of  decomposed 
nodes:  23.6 


Minimum  decomposition 
a=3 

Average  number  of  decomposed 
nodes:  29.4 
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All  a  50%  decomposition 
Average  number  of  decomposed 
nodes:  20.4 


Work  in  progress:  A  Continuous  Wavelet  transform 


During  this  past  study  of  the  fractional  B-splines,  it  became  clear  that  it  would  also  have  been  useful  to  be 
able  to  perform  a  continuous  transform  with  these  functions.  This  has  already  been  done  in  the  case  of 
classical  B-spline  [7],  However,  we  encountered  problems  in  applying  a  similar  scheme  as  in  [7]  in  the  case 
of  fractional  parameter.  More  specifically,  we  need  to  consider  how  can  we  apply  the  formula  (0.14)  in  the 
case  n2  non-integer: 


m  ^ - V - ' 


(^2+1)  times 


where  sequence  associated  to  the  wavelet  spline  functions. 


It  is  more  computationally  practical  to  work  in  Fourier  domain.  Thus,  our  design  approach  to  compute  the 
transform  for  integer  scales  is  in  Fourier.  We  have  written  an  m-scale  relation, 


and  we  can  compute 

P{mco) 

~W) 


0(x/m)  =  Jk)^(x-k) 

keZ 


(  m-X 

e 

it=0 


m 

V  J 


(  m-1 


\a+\ 


=>H,Jen  = 


-yim 


-jko) 


m 

V  J 


Now  that  we  have  the  filter  H,  we  would  like  to  obtain  the  filter  G  that  corresponds  to  the  orthogonal 
wavelet  filter.  Our  approach  to  this  problem  is  illustrated  in  Figure  22  where  hi  and  gi  denote  the  new 
filters. 
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Figure  22.  Filter  bank  for  implementing  a  fractional  B-spline  wavelet  expansion  via  a  2D  CWT. 
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3.  KEY  RESEARH  ACCOMPLISHMENTS 


(1)  During  this  study  we  first  evaluated  exiting  wavelet  functions  that  best  matched  the  smoothness  of  a 
mass  represented  within  a  course  level  of  an  expansion.  We  then  implemented  under  the  Matlab 
programming  environment  an  algorithm  of  the  Continuous  Discrete /racfiona/  spline  wavelet  transform  in 
two  dimensions.  The  algorithm  was  implemented  without  down  sampling,  using  the  “a  trou  Algorithm” 
algorithm  to  provide  redundancy  and  translation  invariance. 

(2)  We  modified  existing  functions  within  the  Wavelab  library  in  order  to  adapt  them  to  the  computation  of 
a  wavelet  packet  representation  and  studied  best  basis  search  strategies  using  the  fractional  spline  wavelet 
transform.  We  have  also  developed  in  interactive  a  tool  for  processing  selected  regions  of  interest, 

(3)  We  ported  this  algorithm  using  the  Matlab  C++  compiler  to  generate  machine  code  that  will  be 
integrated  into  our  PC  based  Mammography  Computer  Aided  Diagnosis  (CAD)  workstation  for  real-time 
interactive  screening  of  mammography  cases  as  described  in  the  Statement  of  Work  for  the  third  year. 
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4.  REPORTABLE  OUTCOMES 


The  exciting  results  of  this  study  will  be  presented  at  the  23rd  Annual  International  Conference  of  the  IEEE 
Engineering  in  Medicine  and  Biology  Society  on  October  25,  2001,  in  a  paper  entitled 
“DETECTION  OF  MASSES  IN  MAMMOGRAPHY  THROUGHREDUNDANT  EXPANSIONS  OF 
SCALE,”  V.  Laffont ,  F.  Durupt,  M.  A.  Birgen,  S.  Bauduin,  and  A.  F.  Laine. 


5.  CONCLUSIONS 


This  phase  of  our  study  showed  that  fractional  splines  functions  are  a  powerful  tool  for  the  representation 
of  masses  in  mammograms  and  are  well  matched  to  the  mass  detection  problem.  The  preliminary  results 
presented  in  this  report  are  promising  but  additional  development  is  needed  during  Task  3.  We  showed  that 
by  using  a  continuous-varying  parameter  for  order,  that  continuous  analysis  outperforms  standard  dyadic 
expansions.  We  showed  that  accurate  approximations  of  mass  shapes  could  be  obtained  through  over 
complete  expansions  of  a  fractional  spline  wavelet  transform.  In  the  context  of  addressing  the  problem  of 
finding  the  best  scale  (as  mentioned  in  last  year’s  annual  report)  the  continuous  transform  implementation 
using  wavelet  packet  theory  appears  to  be  a  very  efficient  method  in  finding  the  best  scale.  Thus, 
reconstruction  after  thresholding  of  coefficients  within  the  best  basis  will  simplify  the  detection  task. 

The  next  phase  of  this  project  will  focus  on  completing  Tasks  2  and  3.  We  have  recently  taken  steps  to 
incorporate  the  fractional  B-spline  algorithm  into  our  CAD  workstation  to  detect  masses  using  these 
representations.  We  will  describe  the  performance  of  this  approach  in  our  next  report. 
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Matlab  files 


Bspline  directory 


function  w=CDWT2D£rac (alpha, type, scale, im) 

% 

%  CWT2Dfracwds  —  2D  Continuous  Discrete  Wavelet  Transform 
%  using  Fractional  B~Splines 


% 

% 

% 

% 

% 

% 

% 

% 

% 


Usage 

w=CWT2Dfracwds (alpha, type, scale, im) 

Inputs 

alpha  spline  parameter  alpha>-0.5 

type  see  FFTfractsplinef ilters 

scale  scale  at  which  the  transform  will  be  computed 

im  image  2D 

Outputs 

w  coefficients  of  the  transform 


% 


%  See  Also: 

%  FFTfractsplinef ilters ,  simu_alpha_scale 


S=size ( im) ; 
n-S(l)  ; 
nc=S  (2)  ; 


if  n-“=nc 

disp(' ERROR  :  image  must  be  square!') 
break 

end 

Jl=log2 (n) ; 

J2=ceil (Jl) ; 
if  J1-=J2 

disp(' ERROR  :  the  size  of  the  input  image  must  be  a  power  of  two!') 
break 

end 

imfft=fft2 (im) ; 

m=n; 

w=  [  ]  ; 

f=[l:m]  ; 

q=l; 

ag=imfft; 
for  k=l: scale 

[Hl,H2] =FFTfractsplinef ilters (m, alpha, tYpe,q) ; 

G=H1(1, :) ; 

G=l/sqrt (2) .*G; 

H=H1(2,:); 

H=l/sqrt(2) .*H; 
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% computation  on  the  rows 
for  i=l:m 

%low~'pass  filtering 
agl (i, : )=ag{i, : ) .*G; 
%high-pass  filtering 
ahl (i, : )=ag(i, : ) .*H; 

end 

%computation  on  the  columns 
for  j=l;m 


a(  : 

j)  = 

agl ( : , j  )  . 

*G. 

dl( 

.1 ) 

=agl ( : , j ) 

.  *H 

d2  ( 

/  j) 

=ahl ( : , j ) 

.  *G 

d3  ( 

/  j) 

=ahl ( : , j ) 

end 

wl=if f t2 (dl ) ; 
w2=ifft2 (d2) ; 
w3=ifft2 (d3) ; 
wl=wl+min  (min  (wl)  )  ; 
w2=w2+min  (min  (w2 )  )  ; 
w3=w3+min  (min  (w3  )  )  ; 

w=  [wl  w2  w3  w]  ; 
ag=a; 

q=q*2; 


end 

approx=if ft2 (a) ; 

appr ox=appr ox +min (min (approx) ) ; 


w= ( [w  approx] ) ; 


* 


[function 

[[FFTanalysisf liters, FFTsynthesisf liters] =FFTfractsplinef liters (M, alpha, 

type,q) 

%  FFTFRACSPLINEFILTERS  Filters  for  the  various  orthogonal  fractional 
%  spline  wavelet  transforms  (orthonormal  or  semi -orthonormal) . 

% 

[FFTanalysisf liters , FFTsynthesisf liters] =FFTfractsplinef liters (M, 
alpha, type) 

%  Computation  of  the  frequency  response  of  the  low-  and  high-pass 

%  filters  that  generate  the  orthonormal  or  semi -orthonormal  (B- 

spline 

%  or  dual)  fractional  splines  of  degree  alpha,  and  of  type  '  +  ' 

(causal) , 

%  (anticausal)  or  (symmetric). 

% 

%  Input : 

%  M  :  size  of  the  input  signal  =  length  of  FFTfilters  =  2^N 

%  alpha>-0.5  :  degree  of  the  fractional  splines 

%  type  :  type  of  B-splines 

%  =  '+ortho'  (causal  orthonormal,  default),  '-ortho' 

(anticausal 

%  orthonormal)  or  '*ortho'  (symmetric  orthonormal); 

%  =  '+bspline'  (causal  B-spline) ,  '-bspline'  (anticausal 

%  B-spline)  or  '*bspline'  (symmetric  B-spline). 

%  =  '+dual'  (causal  dual),  '-dual'  (anticausal  dual)  or 

' *dual ' 

%  (symmetric  dual) .  The  last  option  is  the  flipped  version  of 

the 

%  B-spline  one. 

%  q  :  relative  to  the  scale:  the  filters  change  from  a  scale  to  the 
other 

%  Output : 

%  FFTanalysisf liters  =  [lowpassf liter ;highpassf liter]  :  FFT 

filter  arrays 

%  FFTsynthesisf liters  =  [lowpassf liter ;highpassf liter ]  ;  FFT 

filter  arrays 
% 

%  See  also  FFTWAVELETANALYSIS ,  FFTWAVELETSYNTHESIS . 

%  Uses  FRACTSPLINEAUTOCORR 

% 

%  Author:  Thierry  Blu,  October  1999 

%  Biomedical  Imaging  Group,  EPFL,  Lausanne,  Switzerland. 

%  This  software  is  downloadable  at  http://bigwww.epfl.ch/ 

% 

%  References: 

%  [1]  M.  Unser  and  T.  Blu,  "Fractional  splines  and  wavelets," 

%  SIAM  Review,  Vol .  42,  No.  1,  pp.  43--67,  January  2000. 

%  [2]  M.  Unser  and  T.  Blu,  "Construction  of  fractional  spline 

wavelet  bases, " 

%  Proc.  SPIE,  Wavelet  Applications  in  Signal  and  Image  Processing 

VII, 

%  Denver,  CO,  USA,  19-23  July,  1999,  vol.  3813,  pp.  422-431. 

%  [3]  T.  Blu  and  M.  Unser,  "The  fractional  spline  wavelet 

transform:  definition  and 

%  implementation,"  Proc.  IEEE  International  Conference  on 

Acoustics,  Speech,  and 
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%  signal  Processing  (ICASSP' 2000) ,  Istanbul,  Turkey,  5-9  June  2000, 
vol.  I,  pp.  512-515  . 

if  alpha<=-0.5 

disp('The  autocorrelation  of  the  fractional  splines  exists  only 

') 

disp('for  degrees  strictly  larger  than  -0.5!') 

FFTanalysisf ilters= [ ] ; 

FFTsynthesisf ilters= [ ] ; 
return 

end 

if  M-=2'^round(log(M) /log(2)  ) 
disp('  ') 

disp('The  size  of  the  FFT  must  be  a  power  of  two!') 
disp ( '  ' ) 

FFTanalysisf ilters= [ ] ; 

FFTsynthesisf ilters= [ ] ; 
return 

end 

nu=0:l/M: (1-1/M) ; 


if  nargin<3 

type= ' +ortho ' ; 

end 

if  length (type) <=1 

error ( [ ' ' ' '  type  ' ' ' '  '  is  an  unknown  filter  type! ' ] ) 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%NEW%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

A=fractsplineautocorr (alpha, nu*q) ; 

A2=fractsplineautocorr (alpha, nu* (2*q) ) ; 

%A2=[A  A]  ;  ' 

%A2=A2(l:2;length(A2)  )  ;  %A2(z)  =A(z'^2) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
nu=nu*q;%%%%%%%%%%NEW%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  type(2)=='o' | type (2 ) == ' 0 ' 

%  orthonormal  spline  filters 
if  type(l)=='*' 

lowa=sqrt (2 ) *abs ( (1+exp ( - 
2*i*pi*nu) ) /2 ) . ^ (alpha+1) . *sqrt (A. /A2) ; 
else 

if  type(l)=='+' 

lowa=sqrt (2) * ( (1+exp (- 
2*i*pi*nu) ) /2) . ^ (alpha+1) . *sqrt (A. /A2) ; 
else 

if  type(l)=='-' 


lowa=sqrt  (2)  *  (  (1+exp  (2*i*pi*nu)  )  /2)  .'^(alpha+1)  .*sqrt  (A.  /A2)  ; 
else 

error ([''' '  type(l)  ''''  '  is  an  unknown  filter 

prefix! ' ] ) 

end 

end 


end 

%higha=exp (-2*i*pi*nu) .*lowa; 
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%higha=conj ( [higha(M/2+ {l:M/2) )  higha(l:M/2) ] ) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%% 

higha=exp (-2*i*pi*nu) .*lowa; 

higha=conj ( [higha(M/ (2*q)+ {1:M-M/ (2*q) ) )  higha(l;M/ {2*q) )]); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%% 

lows=lowa; 

highs=higha; 

FFTanalysisf ilters= [lowa;higha] ; 

FFTsynthesisf ilters= [lows; highs] ; 

else 

%  semi -orthonormal  spline  filters 
if  type(l)=='*' 

lowa=sqrt (2 ) *abs ( (1+exp {-2*i*pi*nu) ) 12)  (alpha+1) ; 

else 

if  type(l)=='+' 

lowa=sqrt (2) * ( ( 1+exp ( -2*i*pi*nu) ) /2) (alpha+l) ; 

else 

if  type(l)=='-' 

lowa=sqrt (2 ) * ( ( 1+exp (2*i*pi*nu) ) /2 ) . ^ (alpha+1) ; 

else 

error  (  [ '  '  '  '  type  *  '  '  '  '  is  an  unknown  filter 

type ! ' ] ) 

end 

end 

end 

higha=exp  {-2*i*pi*nu)  .  *lowa.  *A; 

higha=conj ( [higha (M/2+ (1 :M/2 ) )  higha (1 :M/2 ) ] ) ; 

lows=lowa. *A. /A2; 

highs=higha./ (A2.* [A(M/2+(l:M/2) )  A(l:M/2) ] ) ; 
if  type(2)=='d' | type (2) == 'D' 

%  dual  spline  wavelets 
FFTanalysisf ilters= [Iowa; higha] ; 

FFTsynthesisf ilters= [lows; highs] ; 

else 

%  B-spline  wavelets 
if  type(2)=='b' | type (2) == 'B' 

FFTsynthesisf ilters= [lowa;higha] ; 

FFTanalysisf ilters= [lows ;highs] ; 

else 

error ( [ ' ' ' '  type  ' ' ' '  '  is  an  unknown  filter  type! ' ] ) 

end 

end 

end 
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^function  A=£ractsplineautocorr( alpha, nu) 

%  FRACTSPLINEAUTOCORR  Frequency  domain  computation  of  fractional  spline 
%  autocorrelation.  A=fractsplineautocorr (alpha, nu)  computes  the 

%  frequency  response  of  the  autocorrelation  filter 

A (exp (2 *i*Pi*nu) ) 

%  of  a  fractional  spline  of  degree  alpha.  It  uses  an  acceleration 

%  technique  to  improve  the  convergence  of  the  infinite  sum  by  4 

%  orders . 

% 

%  See  also  FFTSPLINEFILTERS 

% 

%  Author:  Thierry  Blu,  October  1999 

%  Biomedical  Imaging  Group,  EPFL,  Lausanne,  Switzerland. 

%  This  software  is  downloadable  at  http://bigwww.epfl.ch/ 

% 

%  References: 

%  [1]  M.  Unser  and  T.  Blu,  "Fractional  splines  and  wavelets," 

%  SIAM  Review,  Vol.  42,  No,  1,  pp.  43 --67,  January  2000. 

%  [2]  M.  Unser  and  T.  Blu,  "Construction  of  fractional  spline 

wavelet  bases," 

%  Proc.  SPIE,  Wavelet  Applications  in  Signal  and  Image  Processing 

VII, 

%  Denver,  CO,  USA,  19-23  July,  1999,  vol.  3813,  pp.  422-431. 

%  [3]  T,  Blu  and  M.  Unser,  "The  fractional  spline  wavelet 

transform:  definition  and 

%  implementation,"  Proc.  IEEE  International  Conference  on 

Acoustics,  Speech,  and 

%  Signal  Processing  (ICASSP' 2000)  ,  Istanbul,  Turkey,  5-9  June  2000, 

vol.  I,  pp.  512-515  . 

N=100;  %  number  of  terms  of  the  summation  for 

computing 

%  the  autocorrelation  frequency  response 


if  alpha<=-0.5 

disp('The  autocorrelation  of  the  fractional  splines  exists  only 

') 

disp('for  degrees  strictly  larger  than  -0.5!') 

A=  [  ]  ; 
return 

end 

S=zeros (1, length (nu) ) ; 
err= [ ] ; 
err0= [] ; 
for  n=-N:N 

S=S+abs  (sine  (nu+n)  )  .  (2*alpha+2 )  ; 

end 

U=2/  (2*alpha+l)  /N'"  (2*alpha+l)  ; 

U=U-1/N^ (2*alpha+2) ; 

U=U+ (alpha+1) * (l/3+2*nu . *nu) /N^ (2*alpha+3 ) ; 

U=U- (alpha+1) * (2*alpha+3) *nu. *nu/N^ (2*alpha+4) ; 

U=U.  *abs  (sin  (pi*nu)  /pi)  .  (2*alpha+2 )  ; 

A=S+U; 
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[function  A=£ractsplineautocorr (alpha, nu) 

%  FRACTSPLINEAUTOCORR  Frequency  domain  computation  of  fractional  spline 
%  autocorrelation.  A=fractsplineautocorr (alpha, nu)  computes  the 

%  frequency  response  of  the  autocorrelation  filter 

A{exp (2*i*Pi*nu) ) 

%  of  a  fractional  spline  of  degree  alpha.  It  uses  an  acceleration 

%  technique  to  improve  the  convergence  of  the  infinite  sum  by  4 

%  orders . 

% 

%  See  also  FFTSPLINEFILTERS 

% 

%  Author:  Thierry  Blu,  October  1999 

%  Biomedical  Imaging  Group,  EPFL,  Lausanne,  Switzerland. 

%  This  software  is  downloadable  at  http://bigwww.epfl.ch/ 

% 

%  References: 

%  [1]  M.  Unser  and  T.  Blu,  "Fractional  splines  and  wavelets, " 

%  SIAM  Review,  Vol .  42,  No.  1,  pp.  43--67,  January  2000. 

%  [2]  M.  Unser  and  T.  Blu,  "Construction  of  fractional  spline 

wavelet  bases, " 

%  Proc.  SPIE,  Wavelet  Applications  in  Signal  and  Image  Processing 

VII, 

%  Denver,  CO,  USA,  19-23  July,  1999,  vol.  3813,  pp.  422-431. 

%  [3]  T.  Blu  and  M.  Unser,  "The  fractional  spline  wavelet 

transform:  definition  and 

%  implementation,"  Proc.  IEEE  International  Conference  on 

Acoustics,  Speech,  and 

%  Signal  Processing  (ICASSP' 2000)  ,  Istanbul,  Turkey,  5-9  June  2000, 

vol.  I,  pp.  512-515  . 

N=100;  %  number  of  terms  of  the  summation  for 

computing 

%  the  autocorrelation  frequency  response 

if  alpha<=-0.5 

disp('The  autocorrelation  of  the  fractional  splines  exists  only 

') 

disp('for  degrees  strictly  larger  than  -0.5!') 

A=[]  ; 
return 

end 

S=zeros (1, length (nu) ) ; 
err= [ ] ; 
err0= [ ] ; 
for  n=-N:N 

S=S+abs  (sine  (nu+n)  )  .  (2*alpha+2 )  ; 

end 

U=2/ (2*alpha+l) /N'^(2*alpha+1)  ; 

U=U-1/N"^  (2*alpha+2)  ; 

U=U+(alpha+l) * (l/3+2*nu . *nu) /N" (2*alpha+3 ) ; 

U=U- (alpha+1) * (2*alpha+3) *nu. *nu/N^ (2*alpha+4) ; 

U=U.  *abs  (sin  (pi*nu)  /pi)  .  (2*alpha+2 )  ; 

A=S+U; 
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Best  basis  directory 


^function  [bb, value]  =test_best  (iiti/MaxCeep^  alpha) 

%  test_best  --  main  program  for  best  basis  reasearch 
%  Usage 

%  [bb, value]  =  test_best (im,MaxDeep, alpha) 

%  Inputs 

%  im  image  2D 

%  MaxDeep  maximum  depth  of  tree- search 

%  alpha  spline  parameter 

%  Outputs 


%  bb  basis-quadtree  of  best  basis 

%  value  value  of  components  of  best  basis 


% 


%  See  Also 

%  sb_Calc2dStatTree ,  sb_Best2dBasis ,  sb_coef_basis ,  sb_Plot2dPartition 


name_dir= [ ' output_alpha' , num2str (alpha, 3 ) , '_deep' , num2str (MaxDeep, 3 ) ] ; 
mkdi r ( name_di r ) ; 


cd  (name_dir) 
im=f  f  t2  ( im)  ; 


[stats]  =  sb_Calc2dStatTree (im, MaxDeep, alpha, 'Entropy' ,0.2) ; 

[bb  value]  =  sb_Best2dBasis (stats, MaxDeep) ; 
bb(l)=l; 

bb=sb_coef_basis (MaxDeep, bb) ; 
save  best_basis .mat  bb 
bbl=f ind(bb) ; 

disp ([ 'decomposed  nodes  =  ' ,num2str (length (bbl) )] ) 

cd  .  . 


figure 

title ([' alpha=  ' , num2str (alpha, 3 ) , '  MaxDeep=  ', num2str (MaxDeep, 3 )] ) 
ax=axis ;hold  on 
sb_Plot2dPartition(bb^ax, 'y') 
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l^function  [SQTree]  -  sb_Calc2dStatTree  ( img,  D,  alpha,  exit ,  Exit  Par )  ^TFPar, 

%  sb_Calc2dStatTree  --  Put  Wavelet/  Statistics  into  Quad  Tree 
%  Usage 

%  SQTree  =  Calc2dStatTree (img, D, alpha, ent [, Ent Par] ) 


Inputs 

img 

D 

alpha 

ent 


EntPar 

Outputs 

sqtree 


2-d  image,  n  by  n,  n  dyadic 
maximum  depth  of  splitting 
spline  parameter 

type  of  entropy  to  record  in  tree: 


options  are 


'Entropy'  --  Coifman-Wickerhauser  ==>  We  use  this  one 
'Log'  --  sum  log (abs ( th_i ) ) 

'l"^p'  sum  |th_i|^p,  0  <  p  <  2,  p  =  EntPar 

'N(eps)'  --  #>=  eps,  eps  =  EntPar 

'Risk'  --  sum  min(th_i^2,eps'^2)  ,  eps  =  EntPar 

'Sum'  --  sum  th_i 

'SURE'  --  SURE (Thresholding) ,  thresh  =  EntPar 

extra  parameter,  depends  on  type  of  entropy 


quad  tree  filled  entropy  numbers, 

%  SQTree (qnode (d,bx,by) )  contains  entropy  of 

quadlet (d, bx, by, n) 


%  Description 

%  A  packet  table  is  indexed  by  depth,  block  within  depth,  and 

%  coefficient  within  block.  2d  images  require  four 

%  splits  per  parent  block,  we  prefers  to  store  the 
%  entropy  measure  for  the  coefficients  in  a  particular  block 
%  and  the  coefficients  themselves.  sb_Calc2dStatTree 

%  creates  this  statistics  tree  for  the  entropy  measure  chosen 
%  by  the  user. 

%  ! ! !  the  stored  coefficients  are  in  the  Fourier  domain  ! 1 ! 

%  See  Also 

%  sb_Best2dBasis ,  sb_Plot2dBasisTree 

% 


[nl, n2 ] =size (img) ; 
boxcnt  =  1; 

SQTree  =  zeros (1,  2*  (4^D) ) ; 
name= [ ' coef_0_0_0 . dat ' ] ; 
save_volume (img, name, 'float32') ; 
if  nargin  ==  4,  EntPar  =  [ ] ;  end 


for  deep  =  0:D, 


[qmf , FFTsynthesisf ilters] =FFTfractsplinef ilters (nl, alpha, ' +ortho' , 2^^  (de 
ep)  )  ; 

qmf =l/sqrt (2 ) .*qmf; 
for  bx=0: (boxcnt-1)  , 
for  by=0: (boxcnt~l) , 


name= [ ' coef_' , num2str (deep) , , num2str (bx) , , num2str (by) , ' .dat ' ] ; 
A=read_volume (name, [nl  n2  1  1] , 'float32') ; 

object=ifft2 (A) ; 
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ss  =  norm(norm(object) ) ; 
object  =  object  ./(ss); 

SQTree (qnode (deep, [bxby] ) )  = 
sb_CalcEntropy (object , 'Entropy' ,EntPar) ; 

%  We  divide  at  each  scale  the  enropy  by  4  in  order  to 
ensure  the  normalization  of  the  entropy 

SQTree (qnode (deep, [bx  by]))  =  SQTree (qnode (deep, [bx 

by] ) ) /4^ (deep) ; 

disp ( [ ' d  bx  by  = 

' ,  num2str  (deep)  ,  n\am2str  (bx)  ,num2str  (by)  ]  )  ; 

if  deep  <  D,  %  prepare  for  finer  scales 

[ql,q2,q3,q4]  =  sb_DovmQuad (A, qmf ) ; 


namel=[ 'coef_' , num2str (deep+1 ) , , num2str (2*bx) , , num2str (2*by) , ' .da 
t']; 

name2= [ 'coef_' , num2str (deep+1 ) , , num2str (2*bx+l) , , num2str (2 *by) , ' . 
dat '  ]  ; 

name3= [ ' coef_' , num2str (deep+1 ),'_', num2str (2*bx) , , num2str (2*by+l ) , ' . 
dat '  ]  ; 

name4= [ ' coef_' , num2str (deep+1) , , num2str (2*bx+l ),'_', num2str (2*by+l) , 

' . dat ' ] ; 

save_volume(ql,namel, 'float32') ; 
save_volume (q2 , name2 , ' float 3 2 ' ) ; 
save_volume (q3 , name3 , ' float 3 2 ' ) ; 
save_volume (q4 , name4 , ' f loat32 ' ) ; 
clear  ql  q2  q3  q4 
end 
end 

end 

boxcnt  =  boxcnt*2; 
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[function  Ent  =  sb_CaicEntropy(  object  ^Entropy,  par) 

%  CalcEntropy  --  Calculate  entropy  number  of  array 
%  Usage 

%  Ent  =  CalcEntropy (object , ent [, par] ) 

%  Inputs 

%  object  1-d  or  2-d  object 

%  ent  type  of  entropy  to  use:  options  are 

%  'Entropy'  —  Coifman-Wickerhauser 

%  'Log'  —  sum  log(lth_i|) 

%  'l^p'  —  sum  |th_i|^p,  0  <  p  <  2,  p  =  par 

%  'N{eps) '  —  #>=  eps,  eps  =  par 

%  'Risk'  —  sum  min{th_i^2 ,  eps""2)  ,  eps=par 

%  'Sum'  —  sum  th_i 

%  'SURE'  —  SURE (Thresholding) ,  thresh  =  par 

%  par  extra  parameter,  depends  on  type  of  entropy 

%  Outputs 

%  Ent  Entropy  of  object 

% 

%  Description 

%  It  is  traditional  to  normalize  the  object  to  norm  1  before 
%  calling  this  routine.  This  routine  does  NO  pre-scaling. 

%  ' 

%  See  Also 

%  sb_Calc2dStatTree 

% 


if  strcmp (Entropy, 'Entropy' ) , 
p  =  (object ( : )  ) . ^2 ; 

Ent  =  -  sum (sum (p  .*  log(p))); 
elseif  strcmp (Entropy, 'Log' ) 
p  =  abs (object (: )  )  ; 

Ent  =  sum(log (eps+p) ) ; 
elseif  strcmp (Entropy, 'l^p' ) , 

%  par  =  p  =  exponent 
p  =  abs (object (: )  )  ; 

Ent  =  sum(p  .'^par); 
elseif  strcmp (Entropy, 'N (eps) ') , 

%  par  =  eps 
p  =  abs (object (:))  ; 

Ent  =  sum(  p  >  par  ) ; 
elseif  strcmp (Entropy, 'Risk')  , 
p  =  object  (:)  .'^2; 

Ent  =  sum(min(p,par'^2)  )  ; 
elseif  strcmp (Entropy, ' Sum' ) , 
p  =  object ( : ) ; 

Ent  =  sum(p) ; 

elseif  strcmp (Entropy, 'SURE' ) , 

tt  =  par'^2;  %  par  =  threshold 

dim  =  length ( obj ect (:)) ; 

p  =  object ( : ) . ^2 ; 

ngt  =  sum(  p  >  tt) ; 

nit  =  dim  -  ngt; 

sit  =  sum(  p  .*  (  p  <=  tt  ) ) ; 

Ent  =  dim  -  2*nlt  +  tt*ngt  +  sit; 

else 

disp('in  CalcEntropy:  unknown  Entropy  request') 
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fprintf  ( 'Request=«%s»\n'  ,  Entropy) 
disp{'I  only  know  how  to  do:') 

disp('  Entropy,  Log,  l^p,  N(eps),  Risk,  SURE  ') 

end 

clear  object 


[function  bb=sb_coef_basis (D^ basis) 

%  sb_coef_basis  return  the  best  basis  tree 
%  Usage 

%  bb=sb_coef_basis (D, basis) 

%  Inputs 

%  basis  best  basis  reurned  by  sb_Best2DBasis 

%  D  Maximiim  depth  of  decomposition 

%  Outputs 

%  bb  best  basis 

% 

%  Description 

%  sb_Best2dBasis  returns  a  tree  with  1  and  0 

%  but  it  can  have  1  on  a  depth  and  only  0  on 

%  the  depth  before. . , 

%  the  coefficients  that  are  non  terminal  nodes 

%  are  deleted 


b=[]  ; 


for  deep  =  0:D-1, 

for  bx=0:  (2'^deep-l)  , 

for  by=0 : (2 ^deep-1 ) , 


nonterminal  node 


end 


end 

end 


if (basis (qnode (deep, [bx, by] ) )  ==  0)  ,  % 

basis (qnode (deep+1 , [2*bx, 2*by] ) )  =  0; 
basis (qnode (deep+1, [2 *bx+l , 2*by] ) )  =  0; 
basis (qnode (deep+1 , [2*bx, 2*by+l ] ) )  =  0; 
basis (qnode (deep+1 , [2*bx+l, 2*by+l] ) )  =  0; 

end 


bb=basis ; 

find (basis) ; 

for  deep  =  D-1:-1:0, 

for  bx=0: (2^deep-l)  , 

for  by=0:  (2"^deep-l)  , 

if (basis (qnode (deep, [bx, by] ) )  ==  0)  ,  % 

nonterminal  node 


name= [ 'coef_' , num2str (deep+1) , , num2str (2*bx) , , num2str (2 *by) , ' .dat 

']  ; 

delete  (name) 

name=  [  'coef_'  , n\im2str  (deep+1)  ,  ,  num2str  (2*bx+l )  ,  ,  num2str  (2 *by)  , ' .d 

at '  ]  ; 

delete  (name) 
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naine=  [  'coef_' , nTiin2str  (deep+1)  ,  , nuin2str  (2*bx)  ,  , num2str  (2*by+l )  , ' .d 

at '  ]  ; 

delete  (name) 

name=:  [  'coef_' , num2str  (deep+l)  ,  , num2str  (2*bx+l)  ,  , num2str  (2*by+l)  ,  ' 

.dat' ] ; 

delete  (name) 

else 

if  (basis (qnode (deep+1, [2*bx, 2*by] ) )  ==  1) 
name= [ 'coef_' , num2str (deep+1) , , num2str (2*bx) , ,num2str (2*by) , ' .dat 

delete  (name) 

end 

if  (basis (qnode (deep+1, [2*bx+l, 2*by] ))  ==  1) 

name= [ 'coef_' , num2str (deep+1) , , num2str (2*bx+l) , , num2str (2*by) , ' .d 
at']; 

delete  (name) 

end 

if  (basis (qnode (deep+1, [2*bx, 2*by+l] ) )  ==  1) 

name=[ 'coef_' , num2str (deep+1) , , num2str (2*bx) , , num2str (2*by+l) , ' .d 
at '  ]  ; 

delete  (name) 

end 

if  (basis (qnode (deep+1, [2*bx+l,2*by+l] ) )  ==  1) 

name= [ 'coef_' , num2str (deep+1 ) , , num2str (2*bx+l) , , num2str (2*by+l) , ' 

. dat ' ] ; 

delete  (name) 

end 

end 

end 

end 

end 


delete ( ' coef_0_0_0 . dat ' ) 
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